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Abstract. 

The Falkner and Skan boundary layer equation is presented 
and the properties of the solutions being sought are discussed. 
The problem is programmed for с with a digital computer 
by two methods, namely, by the Milne Five-Point Method and the 
Runge-Kutta Method. Several solutions are listed and the two 
methods of solution are compared. The results obtained are 
compared with, and found to confirm, the solutions produced by 


D.R. Hartree. 
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Notation.* 


Es » Es» = 


Р' 5 ми. pin 


+ py ir 


U, Y 


U,V, А” 


С Е 


х,у, 


X,Y Z 


a constant of proportionality 

constants denoting convergence requirements 
a function of the boundary layer coordinate, "1 
derivatives of f with respect to ^i 
step-width in the 7 direction 

a characteristic length of the immersed body 
a constant, related to В 

pressure 

& characteristic Reynolds number, = E 


time 


velocity components in the directions parallel to 
and normal to the boundary 


velocity components in rectangular coordinates 
velocity at the outer limit of the boundary layer 
free-stresm velocity 


length coordinates measured tangential to and 
normal to the boundary 


rectangular coordinates 


body forces per unit mass 


* The notation listed here is not applicable to Appendices B and 0, 
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— Um um mood cm «қта 


d e. L З ~ "о 


— 


m Em 
a constant, — Jat) 
boundary layer thickness 


non-dimensional boundary layer coordinate normal 
to the boundary 


kinematic viscosity 
non-dimensional stream function 


density 
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1, Introduction. 

The Falkner and Skan boundary layer equation, 

ги +гг" да - #12) = 0, (1) 
where f is a function of % , with the boundary conditions, 

f=0, f' =O, at % =0, and f' >1 as Y ro, (2) 
describes fluid motion in the laminar boundary layer formed as a 
result of steady, incompressible, two-dimensional flow past an 
immersed body, such as a sharp-edged wedge, a circular cylinder, or a 
parabolic cylinder. This equation was first derived, in a slightly 
different form, by V.M. Falkner and Miss S.W. Skan [1]!. Equation (1) 
was later solved numerically by D.R. Hartree [2] for several values 
of the parameter 8. 

This report describes two methods by which equation (1) can be 
solved with a digital computer. Solutions have been produced with 
a CDC 1604 Digital Computer by both methods illustrated and a 
comparison of the two methods is given. Solutions for several values 
of В , giving 2(7), Р! (07), #"(71), and £'"(7), are listed in 


Appendix D. 


"Numbers in square brackets refer to the Bibliography. 
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2. Тһе Falkner and Skan Equation. 
Prandtl's boundary layer equation for steady, incompressible 
flow [3 | is, in dimensionless form: 


Зи LE. лае (5) 
Falkner and Skan assumed that 
= ах“ (4) 


4, 


where U, is the fluid velocity at the outer limit of the Вайна 
layer, X is the distance, measured tangential to the boundary, from 
the forward stagnation point, and A is a constant of proportionality. 
It has been found [4],[5] that equation (4) closely describes the 
flow over the forward portion of & sharp-edged wedge of included angle 


TA , where the relationship between 3 and mv is given by 


(5) 





= e у 
and over circular and parabolic cylinders, 


Next define a non-dimensional boundary layer coordinate, 


та. ЭТР Қо», (6) 


— Ч , 
and a stream function, 


=, (7) 
29247 Қ», 





where £ = f(7), a function of ” only. The stream function Y 


satisfies the continuity equation, 


(8) 
ох ду , Ж 
namely, 
ду 234 
= Vz = = • 
LL эр? A (9) 





Now the introduction of equations (4) through (9) into equation 
(3) yields, after ві tion, the Falkner and Skan equation (1). 

For a complete derivation of equations (3) and (1) see 
Appendix A, 

Those portions of the solution of equation (1) which are of 
particular interest are the variation of #'(7) with ” and the value 
of £"(0), that is, the value of the unknown boundary condition at 
N = o. It can be shown that 


РА 
Р! se 9 (10) 
ЕЕ 


that is, the variation of £'(7) gives the velocity profile through 
the boundary layer. 


The value of f"(0) is necessary for the determination of the 


local coefficient of friction at the boundary, since the frictional 


shearing stress at the boundary is proportional to а Я апа 
дм) - (тты, «бә | 11 
oy 720 2 2 X £^(0) ко | 


A particular solution of equation (1) which is of special 
interest is the case of boundary layer separation for which it is 
known that £"(0) = 0. In this case, the value of A becomes the 
unknown condition and the particular value of @ which satisfies the 
equation is an essential part of the solution. 

Under certain conditions, equation (1) is applicable to 
compressible flow кеј 7) . In this case, the solution f(%) is 
required for the solution of the energy equation which is essential in 
the determination of heat transfer coefficients. The case of 


compressible flow is not considered in this report. 





3. The Boundary Conditions. 

In addition to the boundary conditions (2) given with the 
differential equation (1), certain additional requirements were 
imposed upon the solutions illustrated in this report. In Figures 1, 
2, and 3, curves of f!, f", and f"' versus 77 have been plotted for 
various values of the parameter /2 . From an examination of these 
curves, the following conditions will be noted. 

1. The range of f! is O< f'< 1 ; that is, f' is always positive 
and approaches the asymptote f' = 1 from below. 

2. £"*> 0 and f"> Об ав 7?-> œ; that is, f" is always positive 
and approaches the asymptote f" = O from above. 

5. 1"! — 0 from below as 7 — œ; even though f"! may be 
initially positive (as it is for Ø < 0), it must become and remain 
negative and approach the asymptote #'' = O from below. 

In programming the solution of the Falkner and Skan equation (1), 
the conditions noted above were included with the boundary conditions 
which were imposed upon the solution being produced. Failure to 
satisfy any one of these specified conditions caused rejection of the 


solution being generated. 





4. The Methods of Solution. 

Two well-known numerical methods were chosen for the solution of 
equation (1). They are the "Milne Five-Point Method" [8],[ 9] апа 
the "Runge-Kutta Method" [10]. The primary reason for choosing two 
methods to solve this problem, when one only is sufficient, is that 
the absolute accuracy of a solution produced by either method is not 
known and a detailed error analysis is difficult to perform. It is 
felt that if good agreement.can be obtained between solutions by the 
two methods, then the solutions can reasonably be accepted as being 
accurate. 

In any numerical integration scheme, the reliability of the 
results is very much dependent upon the selection of a suitable 
step-width, 

С (12) 
In general, the smaller the step-width, №, the more reliable will be 
the results. However, it must be remembered that one pays a penalty 
in that the number of points which must be evaluated increases as h 
decreases. Even with a high-speed digital computer, this factor must 
be borne in mind. The optimum value of the step-width was taken to be 
the maximum value of h which would give results which were reproduc- 
able to the same (desired) number of significant figures when a 
smaller value of h was used. The solutions listed in this report 
were calculated to eight places of decimals. In the case of the 
Milne Method, the optimum value was h = 0.01, and for the Runge-Kutta 
Method, h = 0,001 was required. This meant that the Milne Method 





produced a solution in approximately one-tenth of the time required 
by the Runge-Kutta Method. 

A detailed description of both methods is given in Appendices 
В апа С. Both the Milne Method and the Runge-Kutta Method were 
programmed in FORTRAN language (111, 12 | and were used with the 
CDC 1604 Digital Computer. The complete programs are described and 


listed also in Appendices B and C. 





5. Convergence Criteria. 

Inasmuch as the solution of the Falkner and Skan equation (1) 
approaches an asymptote as 7% >-оо „ the question of convergence 
criteria naturally arises; that is, what is the cut-off point of the 
solution and what requirements must be fulfilled at the cut-off point? 
Furthermore, in order to preclude the possibility of the computer 
computing indefinitely on any trial solution in the event that the 
trial solution happened to approach a horizontal asymptote, f' < 1, an 
upper limit of the variable 77 must be arbitrarily fixed. 

Figures 1, 2, and 3 show that the various solutions rapidly 
approach their asymptotic value. In the preliminary stage, using 
Hartree's work as a guide,.an upper limit of Y = 10. was chosen and 
this proved to be satisfactory for all solutions produced and did not 
require revision. 

The convergence requirements which were imposed at the cut-off 


point are as follows: 


1. (1-г!)< я, · (15) 
2: f" « E, (14) 
2. - 175 Ez (15) 


The value of Es was arbitrarily chosen to be 5. = 1.0 x 10-8. 
If the exact solution to equation (1) were known, the values of E, 
end Ez corresponding to a given value of E, could be calculated, but, 
since such a solution is not known, the values of E, and Ez must also 


be set arbitrarily. A number of preliminary runs were made in which 


the convergence requirements on f" and f"! were omitted, From the 





Pi £” was 


results so obtained, it was noted that for Е 77 10 
approximately one order of magnitude greater than (1 - f!) and 
- f'" was approximately one order of magnitude greater than f". 
Based upon this observation, it was arbitrarily decided to set 
Бр = 10. Х Бр = 1,0 х 1077, апі Е, = 10. x E, = 1.0 x 10-6. These are 


the convergence criteria which were imposed upon the solutions which 


were produced by the methods of this report. 





6. Results. 

Several solutions which have been produced by the Milne Method 
are tabulated in Appendix D. These solutions are quoted to eight 
places of decimals. A number of the cases given were also solved by 
the Runge-Kutta Method. As mentioned earlier, the values of the step- 
width, h, that were used are 0.01 for the Milne Method and 0.001 for 
the Runge-Kutta Method. In each case compared, the solutions produced 
by the two methods were identical to seven places of decimals, and any 
discrepancy in the eighth-place was usually small. 

Table I on the following page indicates the nature of the 
differences which exist between solutions produced by the two methods. 

Hartree [2] gives solutions of f'(7) for several values of 6 to 
four places of decimals. If the corresponding solutions listed in 
Appendix D are rounded off to four places, they are found to agree 
with Hartree's solutions to within +.0005 to -.0005, with the 
exception of two entries, which are likely to be typographical errors; 
these exceptions are as follows: 


ko я f! (Hartree) f! (Appendix D) 
0.5 1.6 0.8860 0.876008 
2.0 0.8 0.7858 0.7958. . 

These differences of +.0005 to ~.0003 are within Hartree's estimated 

limits of error. 
The Milne Five~Point Method, which makes use of "predictor" and 


"corrector" formulae in sequence, is said to be exact for a seventh 


degree polynomial [9 |; for any other regularly behaved function, the 





Table I. Maximum Differences Between Runge-Kutta Method Solutions 
and Milne Method Solutions. 


+ 


[a | = И 
Parameter | Differences | 
(TO) 8 | «o re еа етер | 
то ЫК кб. 1078 1, х 10-9 1. х 10-8 | 
ӘЗ uo) | 
10 0. x 10791. x 108 0. x 108 1. x 1082. x 1079 
c e E 2. 1. 1. 2. 
~.16 ОЖ 1. 0. I. ШЕ 
EO о. 2. © 4. 5. 
-10 о. 27 m 1. b: 
0.0 | 0. 2; а о. 0. 
ори РО. 1. 0. 9. т. 
0.2 0. 2. 1. 1. 2. 
0.5 | Ој ||. 1. 1. 1. | 
0.4 0. 12 1. 4. д. 
0.5 0. 1. 1. 1. 2. 
0.6 0. 1. |1. |1. 3. 
0.7 0. х 1079 1. т о. |2. 
0.8 о. ES 1. т s. 
0.9 о. 1. o. 1. 5. 
1.0 о. 1. 4. |1. А. 
m 0. |. d. 1 3. | 
1.6 ао № 1: 2 11. | 
2.0 0. је 1. қ: 10. 
2.4 о. 1 P 1. E | 





exact solution lies between the "predicted" value and the "corrected" 
value. Furthermore, for a seventh degree polynomial, seventh order 
differences should vanish, end, for other functions, seventh order 
differences should be small and free from erratic changes in 
magnitude. 

In order to note the behaviour of these items, the Milne Method 
solution was programmed to give, in addition to the final solution, 
the predicted value of f"(7), the difference between the second esti- 
mate of f"(7) and its predicted value, and a table of differences to 
the seventh order for the final solution of f'(4). On the following 
four pages are shown this information ав obtsined for the solutions of 
д = =.195 and (8 = 0.7 . The column headings shown have the 


following meanings: 


Paf , 
рЕ = 2' , 
D2F = f" , 

рат , 


PRED D2F = Predicted value of f" , 
03 = (Second estimate of f") - (Predicted value of f") , 
D7 = Seventh order difference . 

It was found that 03 was, in every case, less than 1. x 10-8. 
Тира of D7 were found to be less than 1. x 10-6 and its behaviour 
was not considered to be erratic. From these indications, it appears 
safe to say that the Milne Method solutions are accurate to at least 


six places of decimals. Furthermore, the close agreement between the 


11 





Milne Method Sample Solution No. 1. 
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FALKNER AND SKAN EQUATION NUMERICAL SOLUTION 
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Table of Differences for Sample Solution No. 1. 
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Milne Method Sample Solution No. 2. 


FALKNER AND SKAN 
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Table of Differences for Sample Solution No. 2. 
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Milne Method solutions and the Runge-Kutta solutions indicates that 


both solutions sre accurate to seven places of decimals. 
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7• Comments оп the Methods of Solution. 

Ав will be noted from the descriptions of the two methods given 
in Appendices B and C, the Runge-Kutta Method is self-starting, that 
is, it requires only the values at the first point to be given, 
whereas the Milne Five-Point Method requires six starting points. 
This property of the Milne Method requires that a rather cumbersome 
set of equations be used to obtain the six starting points before the 
Five-Point Method proper can be used. However, this presents no 
difficulty when using a high-speed digital computer. Аз has been 
pointed out, the Milne Method produces a comparable solution with a 
considerably larger step-width than that which can be used with the 
Runge-Kutta Method (0.01 versus 0.001) and this means that the Milne 
Method produces a solution in approximately one-tenth the time 
required by the Runge-Kutta Method. Starting with an initial estimate 
of £"(0) = 0.0, the Milne Method requires from three to six minutes 
computing time on the CDC 1604 Computer. This difference in running 
time is due primarily to the method which was used to alter the 
initial estimate of £"(0); namely, the value of £"(0) was altered by 
a fixed increment in one particular decimal position at a time, 
starting with the first decimal place and moving one place at a time 
to subsequent decimal places as the solution was built up. 

An attempt was made to let the solutions proceed beyond the 
cut-off point once the solution to the cut-off point had been found; 


that is, to let the solution continue asymptotically to higher values 
of 7%. It was found that this could be accomplished with the Runge- 
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Kutta Method but that it could not be done with the Milne Method. 
Beyond the cut-off point, the Milne Method would proceed asymp- 
totically for only a short distance and then it would depart ee 
rapidly away from the asymptote. An attempt to rectify this property 
of the Milne Method was not pursued. It appears that if it were 
desired to cerry a solution beyond the cut-off point, as might well be 
the case when dealing with the case of compressible fluids, that the 
most advantageous way would be to determine the solution to the cut- 


off point by the Milne Method and then introduce this solution into 


the Runge-Kutta Method. 
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8. Conclusion. 

The results quoted in this report confirm the solutions of the 
Falkner and Skan equation given by Hartree. The methods of solution 
described herein are particularly suitable for use with a modern 
digital computer and the results produced are considered to be 
accurate to seven places of decimals. 

Extension of the methods described to the appropriate case of 


compressible flow should present little difficulty. 
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APPENDIX A 


THE BOUNDARY LAYER EQUATIONS 


The 88851456 description of any general fluid motion requires a 
complex system of equations consisting of: 

1. the Navier-Stokes equations, 

2. the continuity equation, 

5. the equation of state, and 

4. the energy equation. 

In the case of incompressible fluids, the equation of state and 
the energy equation are not required, and the system of equations is 
reduced to the Navier-Stokes equations and the continuity equation. 
Furthermore, these latter equations are in themselves somewhat 
simplified when the fluid is incompressible. The Navier-Stokes 


equations for an incompressible fluid ares 


ди. + — ++ зи + ur да == 
9t и / 22. (Эа. ди ads) 
E е“ 4 22 
X= pax * ES ay > az” 
W uev W _ 
о» жү рг eae) s 
| IV д” 3 А-1 
ES 54 (365 — ET 


t о > > 
5 9 32123 (9% LW, 5) 
2 e e Эа? Эх dy” 22 
and the continuity equation is: 


әу oy әз 
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This report is concerned with the case of two-dimensional, steady, 
incompressible flow. Furthermore, it is assumed that the body forces 


are negligible. The system of equations (A-1) and (A-2) then becomes: 


ди ди _ __ Эр Ju, Yue 
4 ax 7 © y P к + 232) 


9X 
oy 9% — зол | РЕМ 
“бх” аф р ә де ну (би э) 


(А-5) 


апа є 
А. л” | 
ах 2 x (А-А) 
with the boundary conditions: | 

m=O, V=O, at 4=0 апа (эа, ав ре”. : 

Prandtl further simplified these equations by considering the 
order of magnitude of each term and then neglecting the smaller terms. 
Before proceeding, it is convenient to express the equations (A-5) and 
(A-4) in terms of non-dimensional quantities by referring the linear 
dimensions, X and Ч» to a characteristic length, L, of tho body about 
which the fluid is flowing and by referring the velocities, U, %, and 


U, to the free-stream velocity, U,. The kinematic viscosity, V), is 


І 


related to a characteristic Reynolds number, Re, namely, 
ы "m 
Rem a UE 
and the pressure, ^p, is made dimensionless by referring it to the 
quantity РЧ. The same symbols will be used for the dimensionless 
quantities as for their dimensional counterparts. 


Now, rewriting equations (A-3) and (A-4) in dimensionless form 
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г A a. 


Ab ол... в 


— — шт а а! ай = | 
іна Dg ц, "m — фањ --- > дин — 
a же. 0:5 A a о Ен: 


ar al) = - з а-а аа” а ж - | 





= bw «е: le К 


> Y 
m! а 





gives: 


Es Кес ap 
„9% | u Grace 


ox dy Әх. Re dy 
> (4-5) 
34- Br dnd Е. жу 
and 
au ar 
SY : 54 О | (А-6) 


with the boundary conditions: 
“=0 , V=O , at y=0 and (>, ав Pun S 
In order to compare magnitudes, it is assumed that the boundary 
layer thickness, 4 » ls very small compared to the characteristic 
length, L » 1.90, 
«і, 
and therefore the dimensionless boundary layer thickness, УХ L , which 
will be denoted б, » 18 very small compared to unity, 1.0.6, 
5 << 1. 
Further, it may be shown [ 4] in several exact solutions of the Navier- 
Stokes equations that the boundary layer thickness is proportional to 
the кам» root of the kinematic viscosity. Thus Vr, is of the 
order of М e It is also assumed that U, and X are of the order 1, 
1.е., ~ ~ 1 and x~1. Thus даћу ~ 1 and from the continuity 
equation (A-6) it is seen that "аг 1 but 3476, во that Ww ~ ó, and 
therefore аг, е Furthermore, А 1, 9%” 7, , ауа, 
ил б, ^ ~ 6, > and o" e $, е 
Rewriting equations (A-5) J (1-6) with the orders of magnitude 


written below each term enables a comparison to be made: 
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4 ж oo - ” [4] нина - 7 
көкей ій ж Io — 7 ч— уд 
= е ж 
au нез. ë — ә - -— "— E 





O ca eom е 
=- = uL T. 
. o — е а-а ы ў) а 
^ i и" — ^ 

пе”. Эш ү m | 
о LAC гүү 


am <> = = 





— 4 $ 


Ра х. ач 
ши! 5 5, $^ i т. 
de w __ 1 (з | u : 
as dt s ee + a, (HE + аг к 
1 & 6 1 б, 5 5 
Э. | әу | -6 
577 + Sy О (А-6) 


1 1 


Now, if the terms of order 5, or less are neglected, equations 


(A-5) and (A-6) are reduced to: 


ди au _ _ ер L Yu 
“5575 7-7 әх Ке ay (A-7) 
a Ep (4-8) 
ди сло НИ 


with the boundary conditions: 

6-0, v=0O, at Ех and U7 «, as ээ 

Equation (A-8) indicates that within the boundary layer, the 
pressure gradient normal to the boundary is negligible. Therefore, 


since ?05,70 : эрд а, x end equation (А-7) сап be written: 


ди, du _ _ 42 , ! з= 
U x + 7 = 22 + Re 2j? (A-10) 


Equation (A-10) is the non-dimensional form of the two-dimensional, 
steady, incompressible flow boundary layer equation. In the 
dimensional form it is: 
| з. 
и 4 р л a 1 АР у за. (А-11) 
ax 24 Р 4х a 
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The Falkner and Skan boundary layer equation will now be derived 
іп the non-dimensional form from equation (А-10). Falkner and Skan 


[1] assumed that the velocity, ((, , at the outer limit of the 
boundary layer is proportional to some power, MV, of tho length 


coordinate tangential to the boundary, measured from the forward 


stagnation point, i.e., 
U, = ax (А-12) 


where A is & constant. 
At the outer limit of the boundary layer, “=, and therefore, 
from equation (A-10), 


u __ ds 
! 3x = Ax э (А-15) 


Co is 


Au 
and, noting that = Ple , equation (A-10) becomes: 


ou ou = ыд a 
— „== — : А-14 


Now. define a non-dimensional У layer coordinate, Z , and 


a non-dimensional stream function, » such that: 


— e = аб 2 
77 = н. гэ 2 (А-15) 
- / 2. UK Ч 
ў mr а — Ех * 


where f = f(^?l), & function of 7| only. The stream funotion y 


&nd 
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satisfies the continuity equation, (A-9), namely: 


Эў аж 
: - Y” — — е Аг1 
= 9g у эх (4-17) 


The components of equation (A-14) may now be obtained from 


equations (A-12), (A-15), (A-16), and (A-17), namelys 








k= «, г, (А-18) 
22 „УК + | ф УП | Yu, 
— va = — Р! : (A-19) 
да. -- TI = pr + M = Р! , (А-20) 
ди 2 re ‚ (А-21) 
7 
л 2. 
Be = U, E f" —, end мн 
94 
Месна С | (4-23) 
ах 


where the primes denote differentiation with respect to 72. 
Substitution of these components into equation (A-14) yields, 


after simplifications 


222137 КІ ага + (1-2) «о . (4-24) 


"ut 
DUM tu 


Multiply through by — and define 


/2 = Ее (A-25) 
FA + 


Equation (A-24) then becomes: 
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a es 
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СЕ ò ЖІ 0 





f"! ee - #2) = 0 (A-26) 
with the boundary conditions: 
Ра 0, 7! = 0, at Y = 0, and #!-> 1 аз 7-- 02. 
Equation (A-26) is the Falkner and Skan boundary layer equation 


with which this report is concerned. 
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APPENDIX B 
PROGRAMMING THE MILNE METHOD NUMERICAL SOLUTION 


OF THE FALKNER AND SKAN EQUATION FOR THE DIGITAL COMPUTER 


Ві. The Equation. 

For the purpose of notation in this section, the variables Y and 
f will be denoted by x and y respectively. The Falkner and Skan 
equation then becomes: 

у" +уу" + В (1-у!7) - 0 (8-1) 
where y is a function of x only and the primes denote differentiation 


with respect to x. The boundary conditions to be imposed are: 


y=0, y'= 0, atx=0, (В-2) 
y! — 1 from below as x — œ, and (В-5) 
y" —- O from above, y'! —- O from below, as x >, (8-4) 


В2. Description of the Milne Method. 
The application of the "Milne Five-Point Method" [83,L9J to 
the numerical solution of equation (B-1) is illustrated below. 


Rewrite the differential equation (B-1) as: 


y!" = -y y? - A (1 - yi?) (B-5) 
The initial conditions are: 

Ус = y(0) = 0, (B-6) 

у, = у'(0) = 0, (8-7) 

уЗ = y"(0) = an ASSUMED value, (B-8) 

В = a SPECIFIED value. (B-9) 


From equation (B-5): 
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yo' = y™"(0) = -6, (в-10) 


уг ч -Уу Уу! “ (1 = 28)у; у; = 0. (8-11) 
A suitable stepewidth, 
Ви Ха ТСХ? (8-12) 


must also be selected. The selection of a suitable value of h has 
already been discussed 58 page 5. 
STEP 1.: Estimate values of у!! and yj". 
Milne suggests: 
у! = yg -hyb , and (B-13) 
y, У = y + № ys" : (8-14) 
In order to obtain a better. estimate, the writer preferred to extend 


these equations by means of the Taylor Series expansion to: 


y!| - y)" - h yiV + hóyY - 9 yYl , and (B-15) 
2 6 
п m v h* У h^ vi 
J1 = Ус + ћ yl * 2 Yo + Е Ус • (8-16) 


Equation (B-11) gives the expression for yy. Subsequent 


differentiation yields: 


у} - -у„ уї' - 2(1 - В)уу! у" - (1 – 2/6)у" , and (B-17) 
уҮі = -у YY - (3 - 28)! уст - (4 - 6 8)у" у!" . (8-18) 


Substitution of (8-17) and (B-18) and the initial conditions into 


(B-15) and (B-16) yields, after simplification, 


Із айы и 28-19) - vg (2 - 36 )y} » and (B-19) 
У." = = 8 + ы (28 = yr + Ч, (2 - 589», e (B-20) 


Evaluation of equations (B-19) and (B-20) will constitute Step 1. 
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STEP 2.: Estimate values of У» У} Yi? yi» y 4? and y, in 
that order as follows: 


2 
Ya "Y -a" — үт) 4. — (8-21) 
е 2 
y: а у. + xo + 16y™ + Ср +7 yy р (8-22) 
! ! и п 16y” n ће m В-2 ) 
A УО (8-25 
бошу! В ( " + 16у" b (B-24) 
Уу = Уо + ЭУ Yan ІЗІ) ug ME 
ER grey) rh (B-2 
У-1 YQ 7 Z4 Uy * 1655 +91) +1, ур уаш -22) 
2 
у, + у, + 20у! + 16у! + Туз) +7 he (8-26) 


An indication of the accuracy of these equations may ђе obtained 
from the remainder terms on the right-hand side of (B-25) апа (В-26) 
which are, respectively, Ros = — h? y! and Ro = “Ros » where y" is 


180 


evaluated at some point in the interval between x , end x, . 


-1 1 
Evaluation of the remainder terms will not be carried out in this 
solution. 

STEP 5.: Evaluate у апа y, from the differential equation, 
(B-5). 

STEP 4.: Using the values of y! and ул from Step 3, repeat 
Steps 2 and 3 until no change, to the desired number of significant 
figures, occurs in the calculated values. At this point, the writer 
required agreement to within 1 x 10-8. 


STEP 5.: Estimate values of У2» Yo» and y as follows: 


„Шум -yM + булн) = 2h? yiv , (в-гт) 


= И 
Y, Yo 5 
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y^ = ye + м + hy; + Yo) » and | (В-28) 
ћ 
Уо = Уо + Y + + Y) > (B-29) 


4 
STEP 6.: Evaluate У)! from the differential equation, (8-5). 


STEP 7.: Correct the values of у^» y^» епі у ез follows: 


ур = У, + aye" Л, (8-20) 
ау. + уз + дуу + уд) , and (8-28) 
Уз = у + Ar + hyi + y3). (B-29) 


STEP 8.: Evaluate y>” from the differential equation, (8-5). 

STEP 9.: Using the value of y^ from Step 8, repest Steps 7 
and & until no change, to the desired number of significant figures, 
occurs in the calculated values. Here the writer required agreement 
to within 1 x 1079. 


STEP 10.: Estimate values of JLo» y!» and y , 88 follows: 


yoe w- Sv" тут ум) асуы (8-21) 
yla. y) > у", + hy", з уй) , and ^— (8-32) 
у= У, = — "ау чу) > (8-22) 


ӨТЕР 11.: Evaluate у" from the differential equation, (8-5). 


STEP 12.1 Correct the values of У о» У! о» апа у_ ав follows: 


у", = у" - — + ћуп ут), (B-34) 
ея. е Хэг, + ћу", + уд) » and (8-52) 
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У-2 * Yo T Ду + nte (в-55) 


STEP! 13.1 Evaluate y!» from the differential equation, (B-5). 

STEP 14.: Using the value of yi, from Step 13, repeat Steps 12 
and 13 until no change, to the desired number of significant figures, 
occurs in the calculated values. Here also, the writer required 
agreement to within 1 x 1079, 

The remainder terms on the right-hand side of equations (В-27), 


(B-30), (B-31), and (B-34) are, respectively, Ro7 = + — Ho 


‚ апа К = = Ro • 


1 11 
Вот ав ИН, Ван - “гт 


STEP 15.: Now regard the values caloulated by the previous steps 


as first approximations ќи refine them as follows: 


у", = (у - Ah, Я ya 


15 ° 4 
= вота + 288у!! - 52у!" + Зул"), (В-55) 
пат ср ву em 
а ts + 288у!, - 52у) « Зуу) , апа (8-26) 
2 
у. = (у, - Ву -Е у) 


o 1570 5 “o 

- табу!» + 280y1, - 52y] + 295). (B-37) 
Evaluate y!5 from equation (B-5). 
У = (Ye E x» s uy уі) 


= =2—(- 21у + 94у! + 11у - 11у2) » (B-38) 
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D 
(o2 fo: 2 09h oa, 11h 


h = 
- — 21y", + олдун, + "у, - 11у:) y end (B-39) 


| г 
МБ. + уо) 


= 285007 21у!, + Gaby! + “у: - 1195) : (В-40) 


Evaluate y!4 from equation (B-5). 
e 

— к 19h van 11h v 

"n (+ as 


h 
+ — 
2880 


y + Iba? Ab + gym Е 21у2') , (в-41) 
2 
у = (у + ^ Y p yi") 


h й Ша " 
— у + ТУ + Ау 21y5) » end (в-42) 


2 
. 19h 24 1185. 
Y: (y, * 5 y, * таг ус) 


h 
+ E У + 144у!, + guy, - 21у5) : (B-43) 


1 


Evaluato y from equation (B-5). 


2 
у; = (уя + ЗВ ут - 5 уі”) 


15 
+ wat Sys - 32у"! » 288у!! » 55у") Е (B-44) 
у; = (у, + = Ya M yn) 
+ сре Ву! = Зеу! + 288у1 + 5592) » апд (B-45) 
у = (у, „Ву: - у") 
М -2—(3у' - 52у'!, + 288у) + 55у). (8-46) 
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Evaluate D from the differential'equation, (B-5). 

STEP 16.: Repeat Step 15 until no change, to the desired 
number of significant figures, occurs in the calculated values. At 
this stage, and at all succeeding points, the writer required agree- 
ment to within 1 x 10710, 

The remainder terms on the right-hand side of equations (B-35), 
(B-38), (B-41), and (B-44) are, respectively, В 43% e 


22 7560 
Ron na = —_n/ уіх Ку, = = and Rap = -R а 
58 7560 : 41 28! 34 29 ° 


STEP 17.: Predict Уз, Yz? and Yz as follows: 


y; = уз + = л9у! + Д9бут + Цабу! = олдуш а долур) 
2 > 
+ == — (B-47) 
Уу yt + (Туп, + 32y3 + 12y] + 32y + 7y3) , ana (8-498) 
y, = у, + 2(7у!', + 52у! + 12у! + 32y! + Ту!) . (8-508) 
| та ~ о 1 2 5 | 


ӨТЕР 18.: Evaluate Ys from the differential equation, (B-5). 


STEP 19.: Correct the values of Y: Уз» ала У as follows: 


уЗ = у" + тш + 52у!" + 12у}! + 52у}! + Ту"), (В-51а) 
2h n п п 

y; = у! + 250771; + 52у + 12у + 22у5 + Ту) y and (В-49а) 

yg * Ya * А + 52у} + 12у} + 52уз + Ту!) . (B-50a) 


STEP 20.: Evaluate ys" from equation (B-5). | 
STEP 21.3 Repeat Steps 19 and 20 until no change, to the 
desired number of significant figures, occurs in the calculated 


values. 
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The remainder terms on the right-hand side of equations (B-47) 
and (В-518) are, respectively, 57 = + жеді ћ! уіх апа 
64 


ВН iX . 
R51 “Тоо? 


STEP 22.: Predict the values of Ул? Yp» and y, for n= 4 as 


follows: 


h 
у" = уч с + Aym, - ум + дбутд - Щуш + ллуш.) 


(8-48) 

у! = у! 2h тун + Z2y" 5 + 12у! 2 ж 22у 05 Ty") (В-49) 
п п-4 15 n-4 n-3 n-2 n-1 n 

Уат Уд 9 СУДЫ, > 32у1-5 + 12yl.n + Bayın + Туш) (B-50) 


ЭТЕР 25.: Evaluate y'" from the differential equation, (B-5). 


STEP 24.: Correct the values of Yn’ Yp» апі у Рог п = 4 ав 


follows: 
ур = Урд + ТЕ Уде + 72уд.5 + 12уд-2 + 22Уд 1 + Тур! (В-51) 
yi = Уд + —I А + 22yn-3 + 12yp 5 * 20yp 4 * Typ) (8-49) 
Ya = Yah Eor + 52у! _„ + 12у}_„ + 52у! + Ту) (8-50) 


STEP 25.: Evaluate уге from the differential equation, (В-5). 
STEP 26.: Repeat Steps 24 and 25 until no change, to the desired 
number of significant figures, occurs in the calculated values. 


The remainder terms on the right-hand side of equations (B-48) 


2214 „7 „Ах 
+ h! y 
7560 
= = 1550 n? ујх » where yix is evaluated at some point in the 


and (B-51) are, respectively, Rig = and 


851 
interval between x4, < апі ху Рог Ryg and between x, j and x, for Ёст, 
STEP 27.: Having completed Steps 22 to 26 inclusive for n = 4, 


repeat Steps 22 to 26 inclusive for n » 5, 6, 7, ......, until such 
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time as the solution is cut off. 

Details of the conditions under which a solution will be cut off 
are explained below under "The Program", 

If егт оь is cut off for being unacceptable, revise the 


initial estimate of y"(0) and repeat the entire procedure outlined 


above. 


BX. The Program. 

The procedure outlined above was programmed in Fortran language 
and was used with the CDO 1604 Digital Computer. A block diagram of 
the Milne Method program is given on page 42 and the complete program, 
designated "Program Milne", is given on pages 46 through 58. The 
following brief explanation is given as an aid in interpreting the 
block diagram and the computer program. 

For a given value of the parameter 4, there is one solution of 


the equation, 


И паду“ + 861 = у'2) = 0, (B-1) 
which satisfies the imposed boundary conditions, (B-2), (B-3), and 
(B-4), as described in section 3, page 4. As has been noted, an 
essential part of this solution is the value of y"(0). Conversely, 
for a given value of y"(0), there is one solution and a particular 
-value of Д which meets all the specified conditions. As indicated 
earlier, this latter form of solution is especially of interest for 
the case of boundary layer separation, for which it is known that 
у"(0) = O and the value of В in this case becomes the unknown. 


The problem was programmed to provide either of these forms of 
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solution. That is, if the value of 8 is specified, the appropriate 
solution and the value of y"(0) are found. Alternatively, the value 


об у"(0) may be specified, in which case the appropriate solution and 
the value of 6 are found. 


The program is written to accept a range of values of the step- 
width, h, but is designed to output the solution at intervals of 
Ха Си“ 0.10. This is accomplished by a loop within a loop (the 
"N" loop within the "I" loop) and for this reason, the value of h 
must be chosen such that the value of 0.10/h is an integer. 

It will be noted on the block diagram that there are, in sequence, 
six decisions which test the acceptability of each point calculated. 
The first three are the imposed boundary conditions, and failure to 
fulfil any one of these tests will result in the solution being cut 
off and rejected. Rejection of the solution results in an alteration 
of the estimate of y"(0), (or @ if appropriate), and the commence- 
ment of a new solution. It should be noted that if y' is initially 
positive, the test y'" 1 О is suppressed until such time as y'' first 
becomes negative. 

The last three decisions apply the convergence criteria. Failure 
to meet the first convergence test results in the calculation of the 
next point of the solution, unless an exit is made from the "I" loop 
ht the test I : 101. This latter test checks the current value of x; 
if I = 101, then x = 10. and the solution is rejected, the assumption 
being that all solutions will converge before x = 10. Failure to 


fulfil the last two convergence tests results in rejection of the 





solution. If all three convergence tests are fulfilled, the solution 


is acceptable and complete. 


It will be noted that certain items in the block diagram are 


marked with asterisks. The following notes are applicable. 


ж 


жж 


BLOCK DIAGRAM NOTES: 

No. 1.1 The method which is used to alter the value of y"(0) or 
Ê» namely by definite increments, requires that: 

(1) For a given (fixed) value of ДВ, the initial estimate of 
y"(0) must be less than the final (correct) value of y"(0). 

(ii) For a given (fixed) value of y"(0), the initial estimate of 
D must be greater than the final (correct) value of Д. 

No. 2.3 If у!" is initially positive, this test is suppressed 
until such time as y"' first becomes negative. 

No. 5.: Тһе value of Д is altered if it is desired to find the 


solution for a given (fixed) value of y"(O). 
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Figure 4. 
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То assist the reader in deciphering the program, some of the 
notation used in the program is explained below. 

1. The functions listed at the beginning of the program 
(e.g. Y3DRVF, THREEF, FOURF, etc.) are compact forms of the ЕЕ 
used in the Milne Method. Where possible, these functions make use of 
the similarity of a set of equations to reduce the number of function 
definitions; e.g. the function YCORRF represents the three equations 
(B-49), (B-50), and (B-51). 

2. The variable y is denoted in the program by both the 
variables W and Y, the variable W being used in the "N" loop, the 
variable Y being used in the "I" loop. The following list will 
illustrate this plus the р of subscripting the variables. The 
illustration is based upon a step-width, h = 0.01, i.e., the number 


of steps in the "N" loop is NO = 0.10/h = 10. 
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5. The meanings of some of the additional variables are given 


below. 
рч, DY = у' 
D2W, D2Y ж у" 
D3W, D3Y . 0 уш 
р2үр = y"(0) =  DeY(1) 
глі “УО 
PD2W | = Predicted value of y" 
ч05 ж (Second estimate of y") =- (Predicted value of y") 
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APPENDIX C 
PROGRAMMING THE RUNGE-KUTTA METHOD NUMERICAL SOLUTION 
OF THE FALKNER AND SKAN EQUATION FOR THE DIGITAL COMPUTER 

C1. The Equation. 

For the purpose of notation in this section, the variables N, and 
f will be denoted by x and y respectively. The Falkner and Skan 
equation then becomes: 

y" «y y" « Й(1-у) = 0, (6-1) 
where у is a function of x only and the primes denote differentiation 


with respect to x. The boundary conditions to be imposed are: 


У = О, у! “ 0, гЁ х е О, (0-2) 
у! з 1 from below as x >œ, and (0-5) 
y" > 0 from above, y''' > O from below, as x »—-“0, (0-4) 


02. Description of the Runge-Kutta Method. 
The application of the "Runge-Kutta Method" [10] to the 
numerical solution of equation (C-1) is illustrated below. 


Rewrite the differential equation (C-1) ав: 


y" а – у yh - B(1- yt) = Fly ,у',ул) (2-2) 


The initial conditions are: 


у Ено) =. Os (0-6) 
у, "у'(0) + 0, (0-7) 
у" = у"(0) = an ASSUMED value, and (6-8) 
В „= а SPECIFIED value. 
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аб іі 1 гял O ъъ = YO шш 
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A suitable step-width, 


Вэ Харт Хо» 


must also be selected. 


already been discussed on page 5. 


(C-9) 


The selection of a suitable value of h has 


Now perform, in sequence, for n = O, the calculations indicated 


below. 
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(6-10) 
(6-11) 
(6-12) 
(6-15) 
(0-14) 
(6-15) 
(6-16) 
(6-17) 
(6-18) 
(0-19) 
(C-20) 
(0-21) 


(0-22) 





1 
" n — — 
ye Y x 5 (Р, + er, + > E Fh) (C-23) 
mo oo pP ! n 0-24 
ys TA "ei — 


Having completed, for n = O, the sequence of calculations 
indicated by equations (C-10) to (C-24) inclusive, repeat the entire 
procedure for n= 1, 2, 5, ......, until such time as the solution is 
cut off. Note that equation (C-24) is, in fact, equation (C-10) for 
the next set of calculations. 

If a solution is cut off for being unacceptable, revise the 
initial estimate of y"(0) and repeat the entire procedure outlined 


above. 


05. The Program. 

The procedure outlined above was programmed in Fortran language 
and was used with the CDC 1604 Digital Computer. A block diagram of 
the Runge-Kutta Method program is given on page 63 and the complete 
program, designated "Program Runge", is given on pages 65 through.71. 

The explanation pertaining to "Program Milne” given in Appendix В 
is equally applicable to "Program Runge” and will not be repeated here. 
The reader is referred to page 59 for that discussion. 

For the sake of clarity, the Block Diagram Notes, which apply to 
those items marked with asterisks in the block diagram, are repeated 


below. 
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BLOCK DIAGRAM NOTES: 

No. 1.: The method which is used to alter the value of y"(0) or 
8 » namely by definite increments, requires that: 

(4) For a given (fixed) value of £ , the initial estimate of 
y"(0) must be less than the final (correct) value of y"(0). 

(ii) For a given (fixed) value of y"(0), the initial estimate of 
8 must be greater than the final (correct) value of ё . 

No. 2.: If y!" is initially positive, this test is suppressed 
until such time &s y!" first becomes negative. 

No. 3.: The value of 8 is altered if it is desired to find the 


solution for a given (fixed) value of y"(0). 
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BLOOK DIAGRAM FOR THE RUNGE-KUTTA METHOD NUMERIOAL SOLUTION 


OF THE FALKNER AND SKAN EQUATION 
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To assist the reader in interpreting the program, the meanings 


of some of the variables used are given below. 
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АРРЕКРІХ р 
TABLES OF NUMERICAL SOLUTIONS OF THE 


FALKNER AND SKAN BOUNDARY LAYER EQUATION 


Numerical solutions of the Falkner and Skan boundary layer 
equation, | 
фи „реп + 2(1 ~ 2'2) = 0 
with the boundary conditions 
f= 0, f' = 0,at 7, = 0, end f?! —- 1 ав 97—02, 

are listed in the following tables for various values of the parameter 
В; f is a function of 7) and the primes denote differentiation with 
respect to N. These solutions have been produced by the Milne Method 


of solution. 


Equivalence of the symbols used in the tables: 
ВЕТА = 8 
D2F(0) = £"(0) 
ETA = 7 
Faf 
DF = Ф! 
D2F = f" 
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EQUATION NUMFRICAL SOLUTION 


FALKNER AND SKAN 
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